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ABSTRACT 

We present the latest developments to the phase modulation method for hnding binaries 
among pulsating stars. We demonstrate how the orbital elements of a pulsating binary star 
can be obtained analytically, that is, without converting time delays to radial velocities by 
numerical differentiation. Using the time delays directly offers greater precision, and allows 
the parameters of much smaller orbits to be derived. The method is applied to KIC 9651065, 
KIC 10990452, and KIC 8264492, and a set of the orbital parameters is obtained for each 
system. Radial velocity curves for these stars are deduced from the orbital elements thus ob¬ 
tained. 

Key words: asteroseismology - stars; oscillations - stars; variables - stars; binaries - stars; 
individual (KIC 8264492; KIC 9651065; KIC 10990452). 


1 INTRODUCTION 

Radial velocities are fundamental data of astronomy. Not only in a 
cosmological context, where the recessional and rotational veloc¬ 
ities of galaxies are of interest, but also in stellar astrophysics. A 
time series of radial velocity (RV) data for a binary system allows 
the orbital parameters of that system to be calculated. However, the 
importance of such data, which are meticulous and time-consuming 
to obtain, creates a large gap between demand and supply. 

In the first paper of this series l lMurnhy et alJbof^ . we de¬ 
scribed a method of calculating radial velocity curves using the 
pulsation frequencies of stars as a ‘clock’. For a star in a bi¬ 
nary system, the orbital motion leads to a periodic variation in 
the path length travelled by light emitted from the star and arriv¬ 
ing at Earth. Hence, if the star is pulsating, the observed phase of 
the pulsation varies over the orbit. We called the method ‘PM’ for 
phase modulation. Equivalently, one can study orbital variations in 
the frequency domain, wh ich lead to frequency modulation (EM; 
IShibahashi & Kurt3l2012h . Similar methods of usin g photometry 
to find binary star s have been developed recently bv lKoenl l l2014h 
and lBalonll ^2014) . though the EM and PM methods are the first to 
provide a full orbital solution from photometry alone. Indeed, the 
application of PM to coherent pulsators will produce RV curves for 
hundreds of Kepler stars without the need for ground-based spec¬ 
troscopy, alleviating the bottleneck. 

The crux of the PM method is the conversion of pulsational 
phase modulation into light arrival time delays, for several pulsa¬ 
tion frequencies in the same star. While the phase modulation is 
a frequency-dependent quantity, the time delay depends on the or¬ 
bital properties, only. Hence for all pulsation frequencies, the re¬ 


sponse of the time delays to the binary orbit is the same, which dis¬ 
tinguishes this modulation fr om other astrophysical sources, s uch 
as mode interaction (see, e.g.. lBuchler. Goupil & Hansenlll99'^ . 

Previously, our approach used numerical differentiation of the 
time delays to produce a radial velocity curve, from which the final 
orbital solution was determined. The RV curves thus obtained were 
sometimes unrealistic due to scatter in the time delays. Recognising 
numerical differentiation as the weakness of the method, we have 
now developed a method of deriving the orbital properties from the 
time delays directly, without the need to convert time delays into 
RVs. It is this method that we describe in this paper. The RV curve 
is produced afterward, from the orbital properties, and is no longer 
a necessary step in the analysis. 


2 TIME DELAY ANALYSIS: METHODOLOGY AND 
EXAMPLE I: KIC 9651065 

Let us divide the light curve into short segments and measure 
the phase of pulsation in each segment. This provides us with 
‘time delays’ (TDs), robs(fn), as observational data, where tn 
(n = 1, 2,...) denotes the time series at which observations are 
available. Fig.[T] shows an example TD diagram (for the case of 
KIC 9651065), where time delays vary periodically with the binary 
orbital period. The TD difference between the maximum and the 
minimum gives the projected size of the orbit in units of light sec¬ 
onds. Deviation from a sinusoid indicates that the orbit is eccentric. 
The TD curve is given a zero point by subtracting the mean of the 
time delays from each observation. The pulsating star is furthest 
from us when the TD curve reaches its maxima, while the star is 
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Figure 1. An example of TD 001 ™ (KIC 9651065) using nine differ- 
ent pulsation modes, including one in the super-Nyquist frequency range 
iMumhv. Shibahashi & Kui^l2013l) . The weighted average is shown as 
solid black squares. 


nearest to us at the minima. The sharp minima and blunt maxima 
in Fig.Q] indicate that periapsis is at the near side of the orbit. The 
asymmetry of the TD curve, showing fast rise and slow fall, reveals 
that the star passes the periapsis after reaching the nearest point to 
us. In this way, TD curves provide us with information about the 
orbit. 

Theoretically, time delay is expressed as a function rth(f) of 
time t and the orbital elements: (i) the orbital period Porb, or equiv¬ 
alently, the orbital frequency VoYh := 1/ Porb, or the orbital angular 
frequency D := 2nVoy:h, (ii) the projected semi-major axis ai sin i, 
(iii) the eccentricity e, (iv) the angle between the nodal point and 
the periapsis toQ and (v) the time of periapsis passage tp- Hence, 
these orbital elements can be determined from the observed TD as 
a set of parameters giving the best fitting rth it). 


2.1 Least squares method 

The best fitting parameters can be determined by searching for the 
minimum of the sum of square residuals 

X^{x, ^) := X] -y - Toba(fn) “ A]^ , (1) 

n 

where (t„ denotes the observational error in measurement of 
'robs(frt)- Here the parameter dependence of Tth(f) is explicitly ex¬ 
pressed with the second argument x, which denotes the orbital el¬ 
ements as a vector, and a parameter A is introduced to compensate 
for the freedom of robs(fra) = 0 (i.e. the arbitrary vertical zero- 
point). Hence the parameters x and A satisfying /dx = 0 and 
dx^/d\ = 0 are to be found, that is, 

-y [rthitn, x) - robs(fn) - A] = Q (2) 

^^ OX 

n ^ 

and 

A=(V-! 5 -) V-ijlrthCfn,®) - robs(fn)]. (3) 


^ We have chosen to represent this angle with zu, rather than with uj, be¬ 
cause of the common use of uj to represent angular oscillation frequencies 
in asteroseismology. 



Figure 2. Top: Schematic side view of the orbital plane, seen from a 
faraway point along the intersection of the orbital plane and the celestial 
sphere, NFN', where the points N and TV' are the nodal points, respec¬ 
tively, and the point F is the centre of gravity of the binary system; that is, 
a focus of the orbital ellipses. The orbital plane is inclined to the celestial 
sphere by the angle i, which ranges from 0 to tt. In the case of 0 ^ T < 7t/2, 
the orbital motion is in the direction of increasing position angle of the star, 
while in the case of n/2 < z ^ tt, the motion is the opposite. The 2 -axis is 
the line-of-sight toward us, and 2 = 0 is the plane tangential to the celestial 
sphere. Bottom: Schematic top view of the orbital plane along the normal 
to that plane. The periapsis of the elliptical orbit is P. The angle measured 
from the nodal point TV, where the motion of the star is directed toward us, 
to the periapsis in the direction of the orbital motion of the star is denoted as 
zu. The star is located, at this moment, at S on the orbital ellipse, for which 
the focus is F. The semi-major axis is a and the eccentricity is e. Then OF 
is ae. The distance between the focus, F, and the star, S, is r. The angle 
PFS is ‘the true anomaly’, /, measured from the periapsis to the stai* at 
the moment in the direction of the orbital motion of the star. ‘The eccentric 
anomaly’, u, also measured in the direction of the orbital motion of the star, 
is defined through the circumscribed circle that is concentric with the orbital 
ellipse. Figure and caption from Shibahashi et al. (2015), this Volume. 


2.2 Time delays as a function of orbital elements 

In order to solve equation we have to derive the ex¬ 
plicit dependence of rth on the orbital elements. The readers 
may c onsult wit h the literature such as Ipreire. Kramer & Lvnel 
(1200ih ('Erratum: l^eire. Kramer & Lvnell2009h . We derive rth fol- 
lowing S hibahashi & Kurtzl (2012h in this subsection. See also 
IShibahashi. Kurtz & Murnh^ (l201.5h . 

Let us dehne a plane that is tangential to the celestial sphere on 
which the barycentre of the binary is located, and let the «-axis that 
is perpendicular to this plane and passing through the barycentre 
of the binary he along the line-of-sight toward us (see Fig.O. The 
orbital plane of the binary motion is assumed to be inclined to the 
celestial sphere by the angle i, which ranges from 0 to 7t. The orbital 
motion of the star is in the direction of increasing position angle, if 
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0 ^ i < 7t/2, and in the direction of decreasing position angle, if 
7t/2 < i ^ 7T. 

Let r be the distance between the centre of gravity and the star 
when its true anomaly is /. The difference in the light arrival time, 
r, compared to the case of a signal arriving from the barycentre of 
the binary system is given by 

r = —r sin(/+ w) sini/c (4) 


where vj is the angle from the nodal point to the periapsis, i is the 
inclination angle, and c is the speed of the light (see FigO. Note 
that r is defined so that it is negative when the star is nearer to us 
than the barycentre0The distance r is expressed with the help of a 
combination of the semi-major axis ai, the eccentricity e, and the 
true anomaly /: 


ai(l - e^) 
1 -I- e cos / 


(5) 


Flence, 


_ ctisini^^ 2^sin/cosrc7-I-cos/sinro 
t(£, X) (1 e j ^ j j. . (o) 


1 -I- e cos / 


The trigonometric functions of / are expressed in terms of a se¬ 
ries expansion of trigonometric functions of the time after the star 
passed the periapsis with Bessel coefficients: 


cos f = —e + 


2(1-e^) 
e 


oo 

Jn{ne) cos ntl{t — tp), 

n=l 


(7) 


Table 1. Observational constraints for KIC 9651065. 


Quantity 

Value 

Units 

"^max 

136 ±27 

s 

"^min 

-211 ±42 

s 

l^orb 

0.003685 ±0.000011 

d-1 

Al 

167.1 ±3.06 

s 

A 2 

35.7 ±3.06 

s 

^('^max) 

0.54 ±0.02 



0.08 ±0.02 




Frequency (d'^) 


Figure 3. Fourier transform of the TD cui've of KIC 9651065 shown in 
Fig.n After identification of Al and A2 in the Fourier transform, their ex¬ 
act values and uncertainties are determined by a non-lineai‘ least-squai'es fit 
to the time-delay cuiwe. 


sin / = Jn'ine) — tp), (8) 

n=l 

where Jn (x) denotes the Bessel function of the first kind of integer 
order n and J„'{x) := dJrt(a:)/da;. Equation with the help of 
equations 0 and 0 gives the time delay rth at time tn for a given 
set of a: = (fl, ai sini/c, e, xu, tp). 

2.3 Simultaneous equations 

Equation 0 forms a set of simultaneous equations for the unknown 
X with the help of equation 0. Let us rewrite symbolically equa¬ 
tion 0 as 

V{x) ■- ^ 


^ That convention is established as follows. When the star lies beyond the 
barycentre, the light arrives later than if the star were at the barycentre: it is 
delayed. When the star is nearer than the barycentre, the time delay is neg¬ 
ative. A negative delay indicates an eaiiy airival time. Since the observed 
luminosity, L, varies as 

L ~ cos(jj{t — d/c), 

where d is the path length travelled by the light on its way to Earth, then 
the phase change, A0, of the stellar oscillations is negative when the time 
delay is increasing. That is, 

T OC —A0. 

The c onvention we hereby establish differs from that in PM I jMurphv et alj 
|20 i 4 equation 3), where the minus sign was not included. We therefore had 
to introduce a minus sign into equations (6) and (7), there, in order to follow 
the convention that radial velocity is positive when the object recedes from 
us. Hence, while the radial velocity curves in that paper have the correct 
orientation, the TD diagrams there are upside-down. Our convention here 
fixes this. 


where 

■— ^ ['rth(fn, x) Xobs^tn) '^ ] ■ (fd) 

This simultaneous equation can be solved by iteration, once we 
have a good initial guess 



Flence we need a means to obtain a good initial guess x^^K 

2.4 Initial guesses 

2.4.1 Input parameters: observational constraints 

An observed TD curve shows, of course, a periodic variation with 
the angular frequency fl. By carrying out the Fourier transform of 
the observed TD curve, we determine Q, accurately. The presence 
of harmonics (2fl, 3D, ... ) indicates that the time delays deviate 
from a pure sinusoid. Flence the angular frequency D is fairly ac¬ 
curately obtained from the Fourier transform of the observed TD 
curve (Fig.0. Let Ai and A 2 be the amplitudes in the frequency 
spectrum corresponding to the angular frequencies D and 2D, re¬ 
spectively. They are also accurately determined, by a simultaneous 
non-linear least-squares fit to the time-delay curve. By folding the 
observational data {Tobs(fn)} with the period 27t/D, we get the 
time delay as a function of orbital phase, </>„ := D(£„ — to)/{2n), 
where to is the time of the first data point. We then know the orbital 
phases at which robs reaches its maximum and minimum. In the 
case of KIC 9651065, shown in Fig.[T] the frequency spectrum is 
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Figure 4. The amplitude ratio between the two components Ai and A 2 of 
KIC 9651065 provides us with an initial guess of e. The thick horizontal red 
line is the measured 2 A 2 /A 1 , and the thin lines above and below it are the 
uncertainties. 


shown in Fig.[^ and the obtained quantities are summarised in Ta- 
ble0 They are used as input parameters from which initial guesses 
for the orbital parameters are deduced. 


2.4.2 Initial guess for e 

The amplitude ratio between the two components Ai and A 2 
provides us with an initial guess of e dShibahashi & Ku^l2012l : 
iMurphv et alj|20f^ : 


J2{2d°'>) 2^2 

Ji(e(o)) “ Ai ’ 


(13) 


where Ji{x) and J 2 {x) denote the first kind of Bessel function, of 
the order of 1 and 2, respectively. This approximation is justified, 
as the vu dependence on the amplitude ratio is weak. In fact, this 
approximation is good for a wide range of vo. Even in the case 
of e ~ 1, the approximation gives = 0.80 (see Fig.|4ll, from 
which the correct value of e is recoverable. In the cas e of e -C 1, the 
LHS of the above equation is further reduced to ~ e jMurphv et alj 

[ 2011 . 


2.4.3 Initial guess for tp 


The largest and the smallest values of robs are well defined and 
easily identified, as are the epochs of these extrema. Therefore the 
extrema and their epochs are useful for providing initial guesses for 
the remaining orbital elements. 

Firstly, let us see when these extrem a occur. With the help 
of the known laws of motion in an ellipse (iBrouwer & Clemencd 
1 19611) . 

d/ ain(l + ecos/) 
dt ^ ’ 

and 


dr ai fie sin/ 


(15) 


where fl denotes the orbital angular frequency, the time variation 


[cos(/ + uj) + e cos U7]. 


of r shown in equation Q is given by 
dr 1 flai sin i 

dt ~ ~cl/T^W 

Hence, when r reaches the extrema 

cos(f + tu) = —e cos zu, 

therefore 

sin(/ + zj) = ±\/l — I 


(16) 

(17) 

(18) 


Since cdr/df = Urad, the extrema of r correspond to the epochs 
of Urad = 0. Geometrically, in Fig.|2 the extrema correspond to the 
tangential points of the ellipse to lines parallel to NN'. Note that the 
nearer side corresponds to negative time delay while the farther side 
corresponds to positive time delay. Hereafter, the orbital elements 
corresponding to the extremum of the nearer side are written with 
a subscript ‘Near’, and those of the farther side are distinguished 
with a subscript ‘Far’. These two points are rotationally symmetric 
with respect to the centre of the ellipse, O. Hence the eccentric 
anomalies of these two points, UNear and upar, are different from 
each other by 7t radians: 


UPar UNear — 7t. 


(19) 


The eccentric anomaly u is written with fl, e and tp as 

u — Q(t — tp) + 2^^ —Jn{ne) sinnQ{t — tp). (20) 
' n 

n = l 

Since the initial guesses for fl and e are already available, the ec¬ 
centric anomaly u in equation d20b is regarded as a function of 
t with a free parameter tp. The epochs of the extrema of T^hs, 
noted as fNear and fpar, respectively, are observationally deter¬ 
mined. Then, by substituting fNear and fpar into equation ( 120b and 
with a constraint given by equation GUI, 


fl(fpar — t-M ear ) 
oo .. 

+2J2-Mne) 

n — l 

X {sinnf2(fpar — tp) 


Sinnfl(fNear “ tp)} — 7t = 0. 


( 21 ) 


This equation should be regarded as an equation with an unknown 
tp. To get a good initial guess for t^^ , we define 

fp := — (tp to). (22) 

27t 

and 


^(0p) :— 27r((()par f/^Near) 
oo 

-f2y -J„(ne) 
n 

n = l 

X {sin 27tn((/par — <t>p) — sin 27tn(/>Near — <?!>p)} — 7r. 


(23) 


We search for zero points of ^(i/p) for a given set of (^Near, <?)par) 
and e = where (()par and 0Near are the orbital phases corre¬ 
sponding to Tmax and Tmin, respectively, that are already measured. 

As in the case shown in Fig.O there are two roots satisfying 

T- (4°)) = 0, (24) 

one corresponding to the case (A) that the pulsating star in question 
passes the periapsis soon after the nearest point to us, and the other 
corresponding to the case (B) that the star passes the apoapsis just 
before the nearest point to us. It is expected then that the sum of w 
derived from these two solutions is 27t, that is, these two solutions 
are explementary angles. 


© 0000 RAS, MNRAS 000, 000-000 




































PM stars II 5 



“iJp 



Figure 7. Two solutions satisfying = D 2 {'cu) = 0. The line of 

sight is assumed to be peipendicular to NN' and the star is viewed from 
the right-hand side. The left panel is the case (A) that the periapsis in the 
near side to us, while the right panel indicates the case (B) that the periapsis 
in the far side from us. It is cleaidy seen that the solution of the case (A) and 
that of the case (B) are explementary angles. 


Figure 5. ^(0p) {red curve) for KIC 9651065. The zero crossings each give 
an initial estimate for T'(0p^^) = 0. Vertical lines at and ‘^n’ 
show the orbital phases of the furthest and the nearest points, corresponding 
to the maximum and the minimum of the time delay, respectively. The light 
cyan bands show the uncertainty ranges of ^p and • 



varpi (rad) 



Figure 6. Discriminants from equations {21\ and i28t for zu of 
KIC9651065. The value of m satisfying both of D\{uj) = 0 and 
D 2 {vj) = 0 can be the solution. The upper panel is the case (A) that the 
periapsis in the near side to us, while the lower panel indicates the case (B) 
that the periapsis in the far side from us. It is clearly seen that the angle zu 
of the case (A) and that of the case (B) are explementaiy angles. 


2.4.4 Initial guess for zj 

Once (f)p is determined, equations 0 and 0 give the true anomaly 
at the nearest point, /Near; 


-f(O) 

cos /ieir = -e 


2 f 1 6^ 1 

H-V Jn{ne) cos 27Tn((/Near “ </p), 

f J 


(25) 

OO 

sin/N°iar = 2-\/l - V Jn{ne) Sin27tn(</)Near - 0p). (26) 

n = l 

Since equations ( 117b and dlSb should be satisfied at the nearest 
point, we define two discriminants 

Dl (vj) := COs(/Near + tu) + 6 COS VJ (27) 


D 2 {w) ~ sin(/Near + Tu) — \/l — 6^ COS^ W (28) 

and search for satisfying both of = 0 and 

D 2 ) = 0 (Fig®- Corresponding to the presence of two pos¬ 
sible solutions of (j>p, there are two solutions of zu, which are ex¬ 
plementary angles (see Fig®. 

It should be noted here that both (j>p and zu are determined as 
functions of e. Their dependence on e for a given set of tf and tn 
is shown in Fig.[8l 


2.4.5 Initial guess for ai sin i 

Once e and zu are determined, the projected semi-major axis, 
ai sini, is determined in units of light travel time, with the help 
of 

'Tmax 'Tmin; by 

aisini (TVnax — Tmin) 2 2 \—1/2 

-= ^^ (1 — e cos zu) . (29) 

c 2 

Note that the two solutions of zu obtained above lead to an identical 
value of ai sini/c. 


2.4.6 TD curve for initial guesses 

Substitution of the initial guesses thus obtained into equation ® 
leads to an initial guess for the TD curve. Among the two possible 
pairs of solutions, one of them generates a reasonable TD curve that 
fits the observations, while the other generates a TD curve that is an 
almost mirror image of the observed TD curve. The value easily 
discriminates between the two values of zu, so this can be auto¬ 
mated. Fig.|^ demonstrates the situation, using the initial guesses 
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Eccentricity 

Figure 8. The dependence of (fip (top panel) and ro (bottom panel) on e for 
a given set of rp and tn of KIC 9651065. In this case, for e < 0.13, there 
is no solution satisfying 'I'(i 7 ip) = 0. It is clearly seen that ro of the case 
(A) and that of the case (B) are explementary angles. 


Table 2. Possible solutions as initial guesses for the orbital parameters of 
KIC 9651065 deduced from the observational constraints listed in Table[T] 
The parameters <j>p and ro given in the first line of each are appropriate to 
be initial guesses, while those in the second line are unsuitable. 


Quantity 

Value 

Units 

^orb 

0.003685 ± 0.000011 


(fli sin 2 )/c 

174 ±25 

s 

e 

0.427 ±0.037 


<l>p 

0.11 ±0.04 

0.51 


zu 

1.90 ±0.23 

4.39 

rad 


for the orbital parameters tabulated in Tablej^ One of the solu¬ 
tions, with periapsis at the far side, fits the data well, while the other 
one having periapsis at the near side has a larger value of /N, 
so the latter is rejected. Of course, the correct solution is consis¬ 
tent with qualitative expectations described in Sect. 1; the periapsis 
of the star is at the near side of the orbit, and the pulsating star 
passes the periapsis after reaching the nearest point to us, that is, 
7t/2 < zu < n. 

Fig.[^ demonstrates how well the TD curve computed from 
the initial guesses reproduces the observed TDs. The orbital phase 



Figure 9. The TD curves for KIC 9651065 constracted with the two sets 
of initial guesses for the orbital parameters. The red curve, generated with 
the parameters in the first line of Table|2] matches the observed ‘time delay’ 
Tobs (violet dots with en'or bar's), wrapped with the orbital period. On the 
other hand, the green curve generated with the pai'ameters in the second 
line of Table[^has a larger value of x^/lV, so it is rejected. The periapsis 
passage 0p was chosen as the orbital phase of zero. The data points Tobs 
were shifted vertically by the amount A so that they match the red curve. 


of zero is chosen so that (j)p = 0. The data points of Tobs were 
vertically shifted by the amount A defined by equation l[3]l, so that 
they match rth. 


2.5 Search for the hest fitting parameters 


Once a set of initial guesses for ai sin i, e, (fip and vj are obtained, 
we may search for the best fitting values of these parameters that 
minimize by iteration. We regard O as a fixed value, because 

the orbital period is already well determined from the Fourier trans¬ 
form of the TD curve. The best fitting values of the orbital param¬ 
eters are summarised in Tablej^ and the TD curve obtained thusly, 
matching best the observed time delay according to the mini¬ 
mization illustrated in Fig.[T^ is shown in Fig.ll II 

The bottom line of Table[^ lists the mass function 
/(mi, m 2 , sin i), defined by 


/(mi,m 2 ,sini) := 


(m 2 sini)® 
(mi -I- m2)2 
{2nfc^ 


G 


2 

"^orb 


ai Sin? 


(30) 


where mi and m 2 denote the masses of the primary (the pulsating 
star in the present case) and the secondary stars, respectively, and 
G is the gravitational constant. The value of /(mi, m 2 , sin i) gives 
a minimum secondary mass of m 2 = 0.82lg Qg, Mq based on an 
assumption of mi = 1.7 M©, so the secondary is probably a main 
sequence G star. 


2.6 Uncertainties 

The uncertainties on the final orbital parameters are the result of 
propagation of the observational uncertainties, which are obtained 
as follows. The uncertainty on the orbital frequency is obtained 
from a non-linear least-squares fit to the TD curve before phase¬ 
folding. 

As for ai sin i/c and e, first, we take a 2-d slice cut of in an 
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Figure 10. ^ jN as a function of (e, ai sini/c) for KIC 9651065. Figure 12. Radial velocity of KIC 9651065. The periapsis passage (pp was 

chosen as the orbital phase of zero. 



Orbital phase 


Figure 11. The best fitting TD curve for KIC 9651065. The periapsis pas¬ 
sage pp was chosen as the orbital phase of zero. 


2.7 Radial velocity 

Since Vmd = cdr/dt, once the orbital parameters are deduced, it 
is straightforward to obtain the radial velocity: 

Oai sini . , . , , 

Wrad =- . [COS(/ + Wj + 6 COS wj . (31) 

vl - e2 

Figm shows the radial velocity curve thus obtained for 
KIC 96 51065. _ 

In lMurphv et alj ( l2014h . we wrote “RV curves derived with the 
PM method could be used as input for codes that model eccentric 
binaries, such as PHOEBE. Given that such codes aim to infer the 
geometry of the orbit, modelling the time delays themselves might 
be preferred over the RV curve, since the former give the binary 
geometry directly and more precisely.^' Our hopes were realized in 
the present work. The radial velocity curve is now provided only as 
a visualisation, —it is not required for the derivation of the orbital 
parameters. 


(e, ai sin i/c)-plane (Fig.llOt. Then by taking a 1-d cut of the plane 
at the values corresponding to the best fitting value of ai sin i/c, we 
get a histogram of along that line. Since the distribution about 
that line is approximately Gaussian, the FWHM of that Gaussian 
gives the uncertainty. The uncertainty thus evaluated on ai sini/c 
for KIC 9651065 is ~ 5 s, and that on e is 0.02. The uncertainties 
on the other parameters are evaluated in the same manner, and they 
are listed in Table[3 


Table 3. The best fitting orbital parameters of KIC 9651065 deduced from 
the observational constraints listed in Table[T] 


Quantity 

Value 

Units 

^OYh 

0.003684 ±0.000011 

d-l 

(ai sin 2 )/c 

183.2 ±5.0 

s 

e 

0.44 ±0.02 



0.14±0.02 


VD 

2.11±0.05 

rad 

/(mi, m 2 , sin 2 ) 

0.0896 ±0.0074 

M0 


3 EXAMPLE 2: KIC 10990452 

Our method is also applicable to pulsators in binary systems with 
short time delays. In this section, we demonstrate KIC 10990452, 
for which the range of variation in time delay is about 1/4 
of the case of KIC 9651065. Fig.[T3 shows the TD curve for 
KIC 10990452. Deviation from a sinusoid indicates that the orbit 
is eccentric as in the case of Example 1: KIC 9651065. However, 
contrary to the case of KIC 9651605, its maxima are sharp and the 
minima are rounded. These facts indicate that periapsis is at the far 
side of the orbit. Fast fall and slow rise reveal that the star passes 
the periapsis after reaching the farthest point from us. Table|4]sum- 
marises the observational constraints for KIC 10990452. 

Substitution of these parameters into equation d23b leads to 
two roots of 41 [pp) = 0, as shown in Fig.[T4l and each solution has 
■uj satisfying D-iiyuj) = D 2 {cu) = 0, as demonstrated in Fig.ll5l 
whose sum is 27t. Initial guesses for the eccentricity, e, and the 
projected semi-major axis, ai sini, are calculated using equations 
( 113b and d29b . 

Substitution of the initial guesses thus obtained into equation 
(|6j leads to an initial guess for the TD curve. Among the two possi¬ 
ble pairs of solutions, the one giving the smaller value of jN fits 
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-400 


55000 55500 56000 56500 

BJD - 2400000 

Figure 13. TD curve of KIC 10990452 obtained from seven different pul¬ 
sation modes. The weighted average is shown as solid black squares. 


Table 4. Observational constraints for KIC 10990452. 


Quantity 

Value 

Units 

'T’max 

63.2 ± 12.6 

s 

^min 

-45.4 ±10.0 

s 

^orb 

0.0081855 ±0.0000142 

d-1 

±1 

49.22 ±1.95 

s 

A 2 

13.14 ±1.32 

s 


0.77 ±0.02 


<P{'^min) 

0.15 ±0.02 



the observations, as shown in Fig.[T^ The other set with the larger 
value of ')d IN is rejected. 

The best fitting parameters are obtained by searching for the 
minimum of IN as a function of (e, ai sin i). Fig.llTlshows a 2D 
colour map of IN. The best fitting parameters are summarised 
in Table|5] and the TD curve generated with these parameters is 
shown in Fig.[T^ Finally, the radial velocity curve is obtained as 
shown in Fig.ll9l 

As seen in the case of KIC 10990452, the present method is 
applicable without any difficulty to pulsators in binary systems 
showing time delay variations of several tens of seconds. Judging 
from the error bars in the observed time delays of Kepler pulsators, 
and from the binaries we have found so far, we are confident that 
the present method is valid for stars showing time delay variations 
exceeding ~ ±20 s. While it may be possible to find binaries with 
even smaller time delay variations, such cases will be close to the 



‘l>p 




Figure 15. Same as Fig.[6] but for KIC 10990452. 



Figure 16. Same as Fig.|9]but for KIC 10990452. The green curve, gen¬ 
erated with one of the solutions of parameters giving the smaller value of 
X^/N, where N denotes the number of data points, fits the data well. On 
the other hand, the red curve generated with the other set of parameters has 
a larger value, so it is rejected. The periapsis passage (pp was chosen as the 
orbital phase of zero. 


Figure 14. Same as Fig.[^ but for KIC 10990452. 
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Figure 17. x'^ as a function of (e, ai sin i/c) for KIC 10990452. Figure 19. Radial velocity of KIC 10990452. The periapsis passage (/)p was 

chosen as the orbital phase of zero. 


Table S. The best fitting orbital parameters of KIC 10990452. 


Quantity 

Value 

Units 

^orb 

0.008190 ±0.000014 

d-l 

(ai sin 2 )/c 

61.3±8.0 

s 

e 

0.55 ±0.03 



0.89 ±0.02 


VD 

5.81 ±0.05 

rad 

/(mi, m 2 , sinz) 

0.01658 ±0.00649 

M0 


noise level of the data and may require external confirmation. The 
noise limit is discussed further in §|^ 



55000 55500 56000 56500 

BJD.2400000 


Figure 20. TD curve of KIC 8264492 obtained from seven different pulsa¬ 
tion modes. The weighted average is shown as solid black squares. 



Orbital phase 

Figure 18. The best fitting TD curve for KIC 10990452. The periapsis pas¬ 
sage (pp was chosen as the orbital phase of zero. 


4 EXAMPLE 3: THE MORE ECCENTRIC CASE OF 
KIC8264492 

Fig-HO] shows the TD curve for another star, KIC 8264492. Devia¬ 
tion from a sinusoid indicates that the orbit is highly eccentric, and 
the number of harmonics to the orbital period visible in Fig. El as 
well as their high amplitudes, confirm this. Let us see if our method 
is valid for such a highly eccentric binary system. Its maxima are 
sharp and the minima are rounded, indicating that periapsis is at 
the far side of the orbit. Fast fall and slow rise reveal that the star 
passes the periapsis after reaching the farthest point from us. 

The orbital frequency, the amplitudes of the highest compo¬ 
nent and the second one, and the orbital phases at the maximum and 
the minimum of the TDs are deduced from the Fourier transform. 
They are summarized in Table|^ Substitution of these parameters 

Table 6. Observational constraints for KIC 8264492. 


Quantity 

Value 

Units 

'^max 

214.6 ±42.9 

s 

’^min 

-132.5 ±26.5 

s 

l^orb 

0.0039408 ±0.0000158 

d-1 

Ai 

159.90 ±3.61 

s 

^2 

50.41 ±3.61 

s 

^('^max) 

0.76 ±0.02 



0.12 ±0.02 
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0.00 0.01 0.02 0.03 0.04 0.05 

Frequency (d'^) 


Figure 21. Fourier transform of the TD curve of KIC 8264492 shown in 
Fig. |2Q| Exact multiples of the orbital frequency (0.00394d“^) are indi¬ 
cated, showing the many harmonics and implying high eccentricity. 



“ijp 


Figure 22. Same as Fig.[^ 'E'(</)p) (red curve), but for KIC 8264492. 




Figure 23. Same as Fig.|6]but for KIC 8264492. The upper panel is the case 
that the periapsis is in the near side to us, while the lower panel indicates 
the case that the periapsis is in the far side from us. 


into equation l l23t enables numerical root-finding of = 0, 

as shown in Fig.[^ One root corresponds to the case (A) that the 
pulsating star in question passes the periapsis soon after the nearest 
point to us, and the other corresponds to the case (B) that the star 
passes the apoapsis just before the furthest point from us. Corre¬ 
sponding to the presence of two possible solutions of there are 
two solutions of ro, which are explementary angles (see Fig. l23b . 

Substitution of the initial guesses thus obtained into equation 
© leads to an initial guess for the TD curve. As in the case of 
KIC 9651605, among the two possible pairs of solutions, one of 
them generates a reasonable TD curve that fits the observations, 
while the other generates a TD curve that is an almost mirror image 
of the observed TD curve, with a larger value of jN (Fig. l24b . 
The correct solution is consistent with qualitative expectations de¬ 
scribed at the beginning of this subsection; the periapsis of the star 
is at the far side of the orbit, and the star passes the periapsis after 
reaching the farthest point from us, that is, 37t/2 < w < 2n. 

Fig.l24l shows the TD curve, computed for the initial guesses, 
plotted with the observed TDs. The orbital phase of zero is chosen 
so that (pp — 0. The data points of robs were vertically shifted by 
the amount A, defined by equation so that they match rth- Un¬ 
like the earlier example of KIC 9651065, there remain systematic 
residuals in the TD curve for KIC 8264492. 

The best fitting parameters are summarised in Tablej?] and 
Fig-H^ and Fig.|^ show the TD curve and the radial velocity 
curve generated with these parameters, respectively. Hence, with 



Figure 24. Same as Fig.|9]but for KIC 8264492. The violet curve, generated 
with the parameters in the first line of Table[6] fits the data well. On the other 
hand, the red curve generated with the other set of parameters has a larger 
value of so it is rejected. The periapsis passage <j>p was chosen as 

the orbital phase of zero. 
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Figure 25. ')p' /N as a function of (e, ai sin i/c) for KIC 8264492. 



Orbital phase 

Figure 26. The best fitting TD curve for KIC 8264492. The periapsis pas¬ 
sage cj>p was chosen as the orbital phase of zero. 


Figure 27. Same as Fig.[T2]but for KIC 8264492. 


5.1 Improvement in technique 

The improvement can be seen in two ways. Firstly, the quality of the 
fit of the theoretical time delay curve to the data can be evaluated 
in terms of the reduced parameter. The former method gave a 
value of 2.21 for KIC 9651065, compared to 1.80 for the analytical 
approach presented in this work. Secondly, one can compare the 
uncertainties in the orbital parameters obtained by each method. 
Table[8] shows that smaller uncertainties, particularly in vo, result 
from fitting the time delays directly, rather than fitting the radial 
velocities obtained by pairwise differences of the time delay data. 

There are also new improvements in the elimination of sys¬ 
tematic errors. Previously, the eccentricity would be underesti¬ 
mated due to the reliance on the approximation in equation ( 113b 
(equation 5 in lMurphv et alj|2014h . This was also strongly subject 
to noise spikes in the Fourier transform of the time delays. Now, 
that approximation is only used as an initial guess, and the search 
for the minimum in the distribution obtains the best-fitting value 
more reliably. 


KIC8264492, we have demonstrated the validity and utility of the 
PM method, even for systems with high eccentricity. 


5 DISCUSSION 

In this work, we have primarily focussed on deriving full or¬ 
bital solutions for highly eccentric binaries. Our first example, 
KIC 9 651065, was also studied in our previous work dMurphv et alj 
l2014h . and so a direct evaluation of the improvement in technique 
is possible. 

Table 7. The best fitting orbital parameters of KIC 8264492 deduced from 
the observational constraints listed in Table|^ 


Quantity 

Value 

Units 

^ovh 

0.003940 ±0.000016 

d-l 

(ai sin 2 )/c 

204.8 ±25.8 

s 

e 

0.67 ±0.04 


0p 

0.80 ±0.02 


zu 

5.28 ±0.05 

rad 

/(mi, m 2 , sin 2 ) 

0.14308 ±0.05410 

M0 


5.2 Factors affecting the minimum measurable time delay 

The detection of the smallest companions, which give rise to the 
smallest time delays, requires a thorough understanding of the 
dominant contributors to the noise and how that noise can be miti¬ 
gated. 

There are many ways that the noise level is affected by the 
properties of the pulsation and/or the sampling. The cadence of 
the observations has little impact on the quoted 20-s limit because 
Kepler observations were mostly made in a single cadence (long- 
cadence at 30 min), though for stars with ample short-cad ence (60- 
s) dat a the phase errors can be reduced by a factor ~ 5 dMurphvI 
|2012|) . The noise can be reduced when the star oscillates in many 
modes, providing they have similar amplitudes to the highest am¬ 
plitude mode. The noise in the weighted average time delays is then 
reduced, though taking the weighted average means that the inclu¬ 
sion of more modes of much lower amplitudes than the strongest 
mode does not help, since phase uncertainties scale inversely with 
amplitude. Also for this reason, we do not consider modes with am¬ 
plitudes below one tenth of that of the strongest mode in each star, 
and high-amplitude pulsators are clearly more favourable. Further¬ 
more, we consider a maximum of nine modes per star due to di¬ 
minishing return in computation time. Finally, it is noteworthy that 
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Table 8. Comparison of the uncertainties in the orbital parameters for KIC 9651065: (i) Those calculated here by fitting the time delay data in this work, vs. 
(ii) those calculated through fitting radial velocities obtained by taking pairwise differences of the time delay data in previous work. 


Quantity 

Units 

Value 

This work 

Value 

Previous work 

^orb 

d-l 

0.003684 ±0.000011 

0.003667 ±0.000016 

(ai sinQ/c 

s 

183.2 ±5.0 

185.0 ±10.0 

e 


0.44 ±0.02 

0.47 ±0.03 

VJ 

rad 

2.11 ±0.05 

2.01 ±0.30 

/(mi,m2,sin2) 

Mo 

0.0896 ±0.0074 

0.0916 ±0.0108 

xVa 


1.80 

2.21 


the smallest detectable time delay variation has no theoretical de¬ 
pendence on the orbital period, providing the orbit is adequately 
sampled. 


6 CONCLUSIONS 

We have developed upon our previous work ^Murphy et alj|201^ . 
where we showed how light arrival time delays can be obtained 
through pulsational phase modulation of binary stars. Formerly, ra¬ 
dial velocities were calculated numerically from the time delays 
and the orbital parameters were obtained from the radial velocity 
curve. Flere, we have shown how the same orbital parameters are 
obtainable directly from the time delays. The radial velocity curve 
is now provided only as a visualisation; it is not a necessary step in 
solving the orbit. 

We will be applying this method to the hundreds of classi¬ 
cal pulsators for which we have measured time delays with Kepler 
data, with the aim of delivering a catalogue of time delay and radial 
velocity curves alongside orbital parameters in the near future. We 
likewise encourage developers and users of binary modelling codes 
to consider taking time delays as inputs. 
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